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Abstract 

We present a new high-order compact scheme for the multi-dimensional Black-Scholes 
model with application to European Put options on a basket of two underlying assets. The 
scheme is second-order accurate in time and fourth-order accurate in space. Numerical ex¬ 
amples confirm that a standard second-order finite difference scheme is significantly outper¬ 
formed. 


1 Introduction 


The multidimensional Black-Scholes model for option pricing (e.g. m) considers n G N >2 
underlying assets Si G [0, oo[ for i = 1,... ,n, where each asset follows a geometric Brownian 
motion, 

dSi{t) = ^iiSi{t)dt + aiSiit)dW^^\t), (1) 


where /ii G R is the drift and (Ji > 0 is the volatility of the asset Si, respectively, for i = 
and dW*'*^(t) denotes a Wiener Process at time t G [0, T] for some T > 0. The 
correlation between the assets is given by dW^^\t)dW^^\t) = pijdt. The Lemma of Ito 
and standard no-arbitrage arguments lead to the following (backward in time) parabolic 
partial differential equation with mixed second-order derivative terms for the option price 
V = ViSi,S 2 ,...,Sn,t) (see, e.g. [S]), 


^ , if- 2c.2dW 
dt + 2 ^"^* * dS^ 

Z=1 '■ 


+ 


dSidSi 

1,3 = 1 
i<j 


^±rS.§.-rV^0, 


with Si > 0, t £ [0, Tl and r > 0 denoting the riskless interest rate. When examining a 
European Put basket option, the final condition is given by 


ViSl,...,Sn,T) =1 




n 

where the asset weights satisfy uii = 1 and additionally oji > 0 for i = 1,..., n if we have 

1=1 

short-selling restrictions. Suitable boundary conditions are discussed later. 

The transformations 


Xi 




T = T-t 


and 


u = 



( 2 ) 
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where 7 > 0 is a constant scaling parameter, yield the (forward in time) parabolic partial 
differential equation 

n 

- 7 ^ ^ Pij 

i<j 

where x € and r G =]0, T]. Under the same transformations the initial condition for a 
European Put basket is given by 

u{xi, ...,Xn,0) =maxf 1 — y^^ujje -r ,0). (4) 

7 i=i 7 

When looking for numerical methods to approximate solutions to problem m, 0 , sub¬ 
ject to suitable boundary conditions, finite difference schemes can be employed, at least for 
space dimensions up to three. Standard discretisations, however, only yield second-order con¬ 
vergence in terms of the spatial discretisation parameter. Alternatively, high-order compact 
schemes can be used which only use points on a compact computational stencil, while having 
fourth-order consistency in space, see for example [nHllElE] and the references therein. A 
drawback is that the derivation of high-order compact schemes (and their numerical stability 
analysis) is algebraically demanding, hence most works in this area restrict themselves to 
the one-dimensional case. An additional complication is present in m in form of the mixed 
second-order derivative terms. 

In a forthcoming paper [3] we derive new high-order compact schemes for a rather general 
class of linear parabolic partial differential equations with mixed second-order derivative terms 
and time- and space-dependent coefficients in arbitrary space dimension n £ N. In the present 
paper we focus on the multi-dimensional Black-Scholes model 0 , 0 . We present a new high- 
order compact scheme which is second-order accurate in time and fourth-order accurate in 
space. To ensure high-order convergence in the presence of the initial condition 0 with low 
regularity we employ the smoothing operators of Kreiss et al. [^. Numerical examples for 
pricing European Put options on a basket of two underlying assets confirm that a standard 
second-order finite difference scheme is significantly outperformed. 


2 

7 ^ “ 


d^u 

dxidx- 


+ 7^] 


r 

(Ti 


^=0 

dxi ’ 


(3) 


2 Discrete two-dimensional Black-Scholes equation 


For the discretisation of 0 with n = 2 we replace the spatial domain by the rectangle 


e = [x 


( 1 ) „(!) 
min’ 


X a^mL] with -00 < X 


(i) ^ „(0 

min ^ "^min 


< 00 for i = 1, 2. On fl, we define the 


grid 

e n I x\l'^ = x^^L + 'ikh, 1 < 4 < iVj,, fc = 1, 2}, (5) 

where h > 0, Nfc £ N and Xmlx = x^/^^ + Nkh for fc = 1, 2. By g 1 ^ we denote the interior of 
(2) 

G\ '. We present the coefficients of a semi-discrete scheme of the form 


il + 1 12 + 1 

t2=*2-l 


at time r for each point x £ Gh for the two-dimensional Black-Scholes equation using n = 2 
in 0. By C 7 ,ij' 2 (a") we denote the approximation of M(x^^\x^y r) after semi-discretisation 
in space with (x^^ x^^^) £ G^^\ 

The general idea underlying the derivation of the high-order compact scheme is to operate 
on the differential equation 0 as an additional relation to obtain finite difference approxi¬ 
mations for high-order derivatives in the truncation error. Inclusion of these expressions in 
a central difference method for equation 0 increases the order of accuracy to fourth order 
while retaining a compact stencil. A detailed derivation of this scheme and a thorough von 
Neumann stability analysis are presented in a forthcoming paper [3]. In the two-dimensional 
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case we obtain the following coefficients 
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as well as 


M7 + l,i2±l =Mii_l,i2=pl = 


Afil±l,i 2 ^ ]^2 ^ 


127 

■,( 2 ) 


,*2 - r. 5 


Afil,i2±l ^ ]^2 ^ 


*72 _ 

2 fJ 2 


127 


Additionally, g{x,T) = 0 for a; € G/, and t € fir- After presenting the high-order compact 
discretisation for the spatial interior we now discuss the boundary conditions. 


3 Discretisation of the boundary conditions 


The first boundary we discuss is S'i = 0 for some i £ {1, 2} at time t £ [0, r[. Once the value 
of the asset is zero, it stays constant over time, see ©• If only one asset reaches its minimum 
value, using 5i = 0 for i £ {1, 2} in the multi-dimensional Black-Scholes equation with n = 2 
leads to the one-dimensional Black-Scholes equation for the asset Sj with j = {1, 2} \ i. One 
can either transform the solution of the one-dimensional Black-Scholes partial differential 
equation using m or derive a fourth-order compact scheme for these boundaries similarly to 
the space interior. If both asset values are minimal, we have 

x^^^ t) = u(x^^^ x^^^ O') 


for r £]0,Tmax] after transforming with @. 

Upper boundaries are boundaries with Si = > 0 with i £ {1,2} at time t £ [0, T[. 

For a sufficiently large we can approximate 


9U(Si,S2,t)| 

dSi 


I 3^=3° 


= 0 , 


( 6 ) 


with Sk € [S'™", for k = {1,2} \ {i}. If only one underlying asset Si reaches its 

maximum value, using m in the two-dimensional Black-Scholes differential equation leads 
to the one-dimensional Black-Scholes differential equation for the underlying asset Sj with 
j = {1, 2} \ {i}. One can either transform the solution of this equation using ([2]) or transform 
the one-dimensional Black-Scholes differential equation using © and derive a fourth-order 
compact scheme for these boundaries. When both underlying assets reach their maximum 
value, we have 


/ m 

U(Xi 


max max 

X2 , 


\ / max max 

r) =u{xi , X 2 


, 0 ) 
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for r €]0,Tniax] after using the transformations ([21l- Since the boundaries behave similar, we 
have 


/ min max \ / min max r»\ / max min \ / max min 

u{Xx ,X2 ,t) =u{x^ ,X2 ,0), u{Xx ,X2 ,t) =u{Xx , X 2 ,0), 

for r €]0 ; "^max] • 


4 Time discretisation 


We use an equidistant time grid of the form t = k At for fc = 0,..., N-r with Nt G N. Using 
a Crank-Nicolson-type time discretisation with step size At leads to 


il+1 12 + 1 

E E 

jl=Il-l 32=*2-l 


^31,02 + 


At 


Kt 


^ 01^32 


+nJ 2 + 

at each point G where only points of the compact stencil are used. By 

denote the approximation of ,Tk)- For the Crank-Nicolson type time 

discretisation this compact scheme has consistency order two in time and four in space. Thus, 
using At G O (h^), leads to fourth-order consistency in terms of the spatial stepsize h > 0. 


= E E 


^ At ^ 

^ 31,32 ^^ 31,32 


5 Numerical experiments 


In this section we present numerical experiments for the Black-Scholes European Puts basket 
option in space dimension n = 2. According to [^, we cannot expect fourth-order convergence 
if the initial condition uo is only in C° (17). In suitable smoothing operators are identified 
in Fourier space. Since the order of consistency of our high-order compact schemes is four, 
we use the smoothing operator $4 (see [5]), given by its Fourier transformation 


l>4(cu) 



This leads to the smoothed initial condition given by 


Uq (xi,X2) 


3h 3h 

-3h -3h 


X, X2 - y) dxdy, 


for any stepsize h > 0, where 'l? 4 (a;) denotes the Fourier inverse of $ 4 ( 0 ;). If uq is smooth 
enough in the integrated region around (* 1 , 12 ) £ 17, we have uo (xi,X 2 ) = uq {xi,X 2 ). Thus 
it is possible to identify the points where smoothing is necessary for a given initial condition. 
This approach reduces the necessary computations significantly. Note that as h —>■ 0, the 
smoothed initial condition uq converges to the original initial condition uo given in Hence 
the approximation of the smoothed problem tends towards the true solution of 

For examining the numerical convergence rate we use the relative f^-error ||[/ref—U ||^2 /||C/ref ||i 2 , 
as well as the Z°°-error ||t/ref — U||ioo, where tUef denotes a reference solution on a fine grid 
and U is the approximation. We determine the numerical convergence order of the schemes 
as the slope of the linear least square fit of the individual error points in the loglog-plots of 
error versus number of discretisation points per spatial direction. We compare the high-order 
compact scheme to a standard second-order scheme, which results from applying the standard 
central difference operators directly in m with n = 2. We use the following parameters, 

(Ti = 0.25, (72 = 0.35, 7 = .25, r = log(1.05), a;i = 0.35 = 1 - L 02 . 

and K = 10. We set the parabolic mesh ratio Ar/h^ = 0.4, but emphasise that neither 
the von Neumann stability analysis presented in [3] nor additional numerical experiments 
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Figure 1: Absolute ^°“-error and relative P-error for two-dimensional Black-Scholes Basket Put 
with smoothed initial condition. 


reveal any restrictions on this relation, indicating unconditional stability of the scheme. We 
use different values pi 2 = —0.8, pi 2 = 0 and pi 2 = 0.8 for the correlation. In Fig. [T] we 
show plots of the Z°°-error and the relative Z^-error. The high-order compact scheme performs 
highly similar for the three different correlation values, the points are almost identical. The 
numerical convergence orders for the high-order compact scheme range between 3.62 and 3.73 
for the Z°°-error, and between 3.87 and 3.94 for the relative Z^-error. The high-order compact 
scheme significantly outperforms the standard second-order discretisation in all cases. 
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